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Blandford's Argument: The Strongest Continuous Gravitational Wave Signal 
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For a uniform population of neutron stars whose spin-down is dominated by the emission of 
gravitational radiation, an old argument of Blandford states that the expected gravitational-wave 
amplitude of the nearest source is independent of the deformation and rotation frequency of the 
objects. Recent work has improved and extended this argument to set upper limits on the expected 
amplitude from neutron stars that also emit electromagnetic radiation. We restate these arguments 
in a more general framework, and simulate the evolution of such a population of stars in the grav- 
itational potential of our galaxy. The simulations allow us to test the assumptions of Blandford's 
argument on a realistic model of our galaxy. We show that the two key assumptions of the argument 
(two dimensionality of the spatial distribution and a steady-state frequency distribution) are in gen- 
eral not fulfilled. The effective scaling dimension D of the spatial distribution of neutron stars is 
significantly larger than two, and for frequencies detectable by terrestrial instruments the frequency 
distribution is not in a steady state unless the ellipticity is unrealistically large. Thus, in the cases of 
most interest, the maximum expected gravitational-wave amplitude does have a strong dependence 
on the deformation and rotation frequency of the population. The results strengthen the previous 
upper limits on the expected gravitational-wave amplitude from neutron stars by a factor of 6 for 
realistic values of ellipticity. 

PACS numbers: 07.05.Tp, 97.60.Gb, 97.60.Jd, 95.55.Ym 
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I. INTRODUCTION 

Continuous emission from spinning neutron stars is a 
promising source of gravitational waves, but so far no 
detections have been reported. This begs the question 
"what is the largest expected ampHtude of the continu- 
ous signal from nonaxisymmetric neutron stars?" There 
might exist a (so far undetected) population of spinning 
neutron stars whose dominant energy loss goes into the 
production of gravitational waves, rather than into elec- 
tromagnetic radiation. These are often called "gravi- 
tars" ; we show later in this paper that gravitars set an 
upper limit on the amplitude of gravitational waves from 
spinning neutron stars that are also emitting electromag- 
netic radiation. 

In 1984 Blandford found a simple analytic relation- 
ship between the expected maximum amplitude of grav- 
itational waves emitted by gravitars and their average 
galactic birthrate. This argument was not pubhshed but 
it is documented by a citation in [T]. The argument was 
recently revised in [2j. This paper revises both the origi- 
nal and the revised Blandford arguments, and shows that 
two key assumptions of these arguments do not hold in a 
realistic galactic model of gravitars. This paper corrects 
the assumptions of the argument, and then investigates 
how the conclusions are affected by this change. 

We stress that while this paper studies the behavior 
of a population of galactic gravitars, it does not make 
a plausibility case for the possible existence of such ob- 
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jects, or study their potential astrophysical implications. 
The study itself is nevertheless interesting, because even 
if gravitars do not exist, they provide a relevant up- 
per bound on the gravitational-wave emission by objects 
(such as rapidly-spinning neutron stars) that do exist. 

The paper is organized as follows. Sec. [ll] reviews 
Blandford's argument and its assumptions. A simple an- 
alytic calculation is used to derive the frequency-space 
distribution of the sources. This allows a sharper state- 
ment of the conclusion and clarifies the dependence upon 
the assumptions. The aim of this paper is to test whether 
these assumptions are fulfilled in a realistic model of our 
galaxy and, if the assumptions do not hold, what the 
consequences are. Sec. Ill describes a numerical sim- 
ulation of the Galaxy, and Sec. |IV| presents results for 
the simulated spatial and frequency distribution of grav- 
itars at the present time. The simulated spatial distri- 
butions do not satisfy the assumptions of Blandford's 
argument: they are not two-dimensional and uniform. 
Sec. |V] uses these simulated distributions to recompute 
the expected maximum gravitational- wave amplitudes 
from gravitars. As shown in 2 , the maximum expected 
gravitational wave amplitudes from gravitars are upper 
limits for the gravitational wave amplitudes from neu- 
tron stars spinning down through combined electromag- 
netic and gravitational-wave emission. Previous work as- 
sumed that all neutron stars are formed with the same 
(high) birth-frequency. Here, the argument is generalized 
to cover a continuous distribution of initial frequencies. 
This is followed by a short conclusion. 

For realistic models of neutron stars, the general upper 
limit on gravitational-wave emission set by considering 
the gravitar case applies for gravitational- wave frequen- 
cies / > 250 Hz. The reader who wants to skip all the 



details and just see the final result is advised to look at 
Fig. [5] which is the main result of the paper. 



II. BLANDFORD'S ARGUMENT: 
ANALYTIC DESCRIPTION 



AN 



A. Frequency Evolution and a First Analysis 

If a rotating neutron star has a nonaxisymmetric 
shape, it will radiate away rotational energy by the emis- 
sion of gravitational waves. It is straightforward to de- 
rive the equations describing the frequency evolution of 
gravitars. Their spin-down due to their nonaxisymmetric 
shape is given by 



i^-'^^^-f. 



(ILl) 



where / is the frequency of the emitted gravitational 
waves, which is twice the spin frequency of the gravitar. 
G is Newton's gravitational constant, / is the momen- 
tum of inertia with respect to the rotational axis, c is the 
speed of light an d e = ^-^^jM jg \}^q ellipticity of the gravi- 
tar. Integrating (ILl I gives the frequency of gravitational 



waves emitted at time t as 



f(t)^{f-^ + p-^eh) 
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, (II.2) 



assuming an initial (birth) gravitational-wave frequency 
/b = / (0) . The constant /? is approximately 



/33 ^5.3 X 



10^**kgm^ 
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kHz. 



(II.3) 



We infer from (II. 2 1 the spin-down timescale 



Tew (e, /) := Pe-T 



4.6 Gyr 



10- 



100 Hz 



(n.4) 

which is the time for a gravitar born at a gravitational 
wave frequency fs^f with ellipticity e to spin down to 
gravitational-wave frequency / via the emission of grav- 
itational waves. These equations allow one to calculate 
the gravitational-wave frequency at the present time for 
any gravitar given its birth frequency, ellipticity, and age. 
This in turn allows one to determine the frequency dis- 
tribution of a population of gravitars. 

The strain amplitude h of gravitational waves emitted 
by a gravitar at distance r to the detector and assuming 
optimal mutual orientation (gravitar sky position given 
by the unit vector orthogonal to the plane of the detector 
arms; gravitar spin axis parallel to this vector) is given 

by 
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Given a model for the spatial distribution of gravi- 
tars, this allows one to determine the distribution of 
gravitational- wave amplitudes. 

Let us begin by giving Blandford's original argument in 
a more complete form than the single-paragraph version 
given in [T]. 

Assume there is a population of galactic gravitars, 
which remain undetected because they do not emit elec- 
tromagnetic waves. Ref. [3] (particularly Sec. 6 and the 
appendix) shows the conditions necessary for a neutron 
star to be a gravitar. Ref. [4 is a simulation of a popula- 
tion of isolated neutron stars accreting matter from the 
interstellar medium and demonstrates that quite a few 
neutron stars may in fact meet these conditions. Taken 
together, these two papers establish a detailed plausibil- 
ity argument for the possible existence of a population of 
gravitars. 

Assume the neutron stars are uniformly distributed in 
a thin two dimensional galactic disk with radius R and as- 
sume that the time between gravitar births in our galaxy 
is constant: Tb ~ SOyrs. Assume that all gravitars are 
born with the same ellipticity e and high birth frequency 
/b. The frequency of each gravitar will then evolve ac- 
cording to (II.2I. 



Consider an interval [/i , /2] of gravitational wave fre- 
quencies with /i, /2 ^ /b- Let us consider wide ranges of 
frequencies corresponding to the broad-band sensitivity 
of modern interferometric detectors. Then, from (II. 2 1 



the time a gravitar will spend in this interval of frequen- 
cies is 



1-12 



/2 



TGw(e,/i) 



(H.6) 



The number of sources in this frequency interval is given 
by A^i2 = ^ if ti2 ^ Tb and depends on the choice of 
the frequency band. Because of the assumed two dimen- 
sionality and uniformity of the spatial distribution the 
average distance r^ to the closest gravitar in this range 
of frequencies can be written as 
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R 



R 



rGw(e,/i) 
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(H.7) 



The formula given in Ref. [Tj agrees with Eq. (II.7I if 
one assumes that the factor in square brackets is of order 
unity, which is the case for the latest generation of broad- 
band interferometric gravitational-wave detectors. 

To estimate the gravitational-wave amplitude of the 
strongest source requires a bit of care. To get the flavor 
of the original argument, consider a one-octave frequency 
band [/, 2/]. The quantity in square brackets in (II. 7 1 is 
15/16 which we approximate as unity. Substituting the 
distance to the closest source (II. 7) into (II. 5 1 and ne- 
glecting the fact that /^ can vary by up to a factor of four 
within the octave, gives the amplitude of the strongest 



source in this one-octave frequency band to be 
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(II.8) 
The Blandford argument is simply the observation that 
this amphtude is (1) independent of the population's de- 
formation e and (2) independent of frequency /. 

Blandford's argument may also be stated in terms of a 
comparison between two different model Galaxies, each 
containing a similar populations of gravitars but each 
having a different (but constant) value of the ellipticity. 



B. Restating Blandford's Argument 

The previous paragraph is a rigorous version of Bland- 
ford's original argument. We now generalize this, build- 
ing on the methods first presented in [2]. 

Let us first define useful quantities to describe a popu- 
lation of gravitars. As before r denotes the distance be- 
tween gravitar and detector, / is the frequency of grav- 
itational waves emitted, and e is the ellipticity of the 
gravitar. 

In this section t measures the age of the gravitars. 
A gravitar with age t = is born at the present time, 
whereas t > for a gravitar born in the past. 

Because optimal mutual orientation of the gravitar's 
spin and the detector's normal axis is assumed, the 
gravitational-wave amplitude depends only on the dis- 
tance r but not on the sky position. Therefore, it is useful 
to define the probability dP^ of finding a gravitar born 
time t ago at the present time in a spherically symmetric 
shell [r, r + dr] around the Sun [11 . The probability can 
be written in terms of a probability density g^ (r, t) as 



dP, = g,{r,t)dr. 



(II.9) 



Moreover, let us define the probability dPf of finding 
gravitars born time t ago with ellipticity e in a present- 
time frequency band [/, / -I- d/] . dPf can be written in 
terms of a probability density gf (e, /, t) as 



dP/ = e/(e,/,t)d/. 



(11.10) 



Note, that the probability densities are normalized by 
J^ dr g, (r, t) = 1 Vt and J^ d/ gf (e, f,t) = l Ve, t. 

For further generalization consider a continuous distri- 
bution of gravitational wave frequencies at birth instead 
of a single, high value. Let dPf^ be the probability of 
the birth frequency being in a band [/o,/o +d/o]. The 
corresponding probability density gf^ (/q) is defined by 



dPf„ = gf„ (/o)d/o, 



(11.11) 



normalized as before. Frequency change by redshift from 
cosmological evolution is neglected since all gravitars con- 
sidered are within our galaxy. 



To link the initial frequency distribution gfg to the 
present-time distribution gf, consider a gravitar with el- 
lipticity e whose current frequency is /, and let /o (e, /, t) 
denote the gravitar's frequency at time t in the past. 
Solving (II. 2 1 for the birth frequency yields 



/o(e,/,i)-(r 



r'e't) 



(11.12) 



The probability density gf can be rewritten in terms 
of the initial frequency distribution gfg by a change 
of variables. The fraction of gravitars in a birth fre- 
quency band [/o, /o + d/o] is the same as the fraction in 
a present time frequency band [/, / -I- d/] , so the identity 
Qf d/ = gfo d/o yields 

gf (£, /, t) d/ = gfg (/o (£, /, t)) ^J^lMA d/ (11.13) 



Qfo{h{£j,t)) 



P 



d/, (11.14) 



from which 



gf{ej,t) = £'/o(/o(£,/,i)) 






(11.15) 



immediately follows. 

To allow for a time-dependent birthrate of galactic 
gravitars, let h (t) be the birthrate as a function of t. 
The number of gravitars dA^ formed during a short time 
interval [t, t + dt] is then dN — n (t) dt. 

The number dN of gravitars in a thin spherical shell 
[r, r + dr] around the position of the Sun, with frequen- 
cies in [/, / + d/], with fixed ellipticity e, formed in a 
time interval [t, t + dt] ago is then given by 

dN = dPr X dPf X dN 

= ft. (r, t) dr X gf (e, /, t)df xfi (t) dt.(II.16) 



From (II. 5 1 it follows that for fixed e and / there is a 



unique, invertible mapping r (h) from the amplitude of 
gravitational waves h to the distance r of the gravitar 
from the Sun. A change of variables from r to h yields 

dr(h) 
dN = g,{r{h) ,t) ^^dhx gf{ej,t)df X h{t)dt 

(11.17) 
and the number M (/i, /2, /imax) of gravitars with a grav- 
itational wave amplitude h ^ /imax in a frequency band 
[/i I /2] and ages t ^ i is given by integration as 

A^(/l,/2,/imax) = 



* r^'^ r°° dr(h) 

dth{t) / df gf{ej,t) / dhg,{r{h),t) ^ > 

Jfl "''imax 



dh 

(11.18) 



Here, the integral over h is performed for a fixed e (as- 
suming the same ellipticity for every gravitar) and fixed 
/ and t. After integrating out the dependence on h the 
follow-up integration over / weights the previous integral 



by the frequency density. The last integration sums the 
distributions from different birth times weighted by the 
galactic neutron star birthrate at that time. 

Before Eq. (11.18 ) is used to rederive Blandford's result, 



let us prove that the frequency distribution from a single 
birth frequency has reached a steady state {dtQf — 0) 
at frequency /, if / (i) < / < /g, where / (i) is given 
by (II.2I. Consider a frequency band [/, / + d/] which 



is wide enough to contain at least one gravitar at all 
times t. li f (t) < / < /b, the constancy of the birthrate 
guarantees that if and only if a gravitar leaves the fre- 
quency band by the lower boundary another gravitar will 
enter the frequency band from higher frequencies. The 
assumption of a steady state is crucial. If the distribu- 
tion has not reached a steady state in a certain frequency 
band there will be no sources in that band and there is 



no contribution to the integral in (11.18 1 



Let us now re-derive Blandford's result by using a den- 
sity function g,- (r, t) — 2r/R^ which describes a popula- 
tion of galactic gravitars uniformly distributed in a flat 
two-dimensional disk with radius R [T2]. Further, as- 
sume a constant birthrate h{t) = — and a single high 
birth frequency /b such that gf^ (/p) = S (/o - /b). In- 



serting r (h) into (11.181 by solving (II. 5) for r we find 



after a slightly technical but straightforward calculation 



M(/i,/2,/l„.ax) 
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The integrations are trivial and yield 
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(11.19) 



(11.20) 



Let us follow [2 and assume a 50% chance of detec- 
tion, corresponding to M = 1/2. One finds a maximum 
gravitational-wave strain [2] 



'''max 
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(IL21) 



This result can also be directly compared with the earlier 
result (11.81 from the cruder analysis, by setting /2 = 2/i 
and setting M = 1. The values of /imax obtained by these 
two different analyses disagree by about 40%, but are 
independent of deformation and frequency /i . 

For a broad band search performed today we assume 
ln(/2//i) ~ 1- Then (11.21 1 gives the largest amplitude 



expected under the assumptions [T3] from galactic grav- 
itars as /imax « 4 X 10"^''. 

For later comparison with the realistic galactic model, 
let us calculate the dimensionless averaged fractional fre- 
quency density gj (e, /) in the population. It is de- 
fined by dP/ — gf {e, /) -4- being the probability to find 
gravitars with a fixed ellipticity e in a frequency band 
[fj + df], 



^/ (£,/) = 



/ 



N, 



tot Jo 



dt'h{t)gf{sj,t), 



(11.22) 



where TVtot = Jndt h is the number of gravitars formed 
during the timespan t. Using the same assumptions as 



for the derivation of (11.20 1 yields 



Qfi^J) = — e 

IB 



*•'-'!-'. 



(11.23) 



For fixed ellipticity the averaged fractional frequency 
density falls off with Z^'* and scales for fixed frequency 
as £~2. 

Let us summarize the assumptions made for this ana- 
lytic calculation of the strongest gravitational-wave sig- 
nal from galactic gravitars. Assume all gravitars are born 
at a single high birth frequency with fixed ellipticity and 
constant birthrate, and reside in a two-dimensional, uni- 
form distribution, i.e. in a thin galactic disk. Assume 
their spin-down is governed by the emission of gravita- 
tional waves as described by Eq. (II.2I. Adopting these 
assumptions and a 50% chance of actual detection, the 
largest amplitude ft-max of gravitational waves emitted by 
galact ic gravitars in a frequency band [/i , /2] is given by 
(|II.21|). Thus, a precise statement about /imax is the fol- 



lowing 

Result: Assume the existence of a popula- 
tion of galactic gravitars with uniform, two- 
dimensional spatial distribution, single, high 
birth frequency, fixed ellipticity e, and con- 
stant birthrate. Choose a frequency band 
[/, sf] with scale s > 1 large enough such 
that there is at all times at least one gravi- 
tar in this band. Then the largest amplitude 
hmax of gravitational waves emitted by galac- 
tic gravitars in this band is independent of f 
and e and depends only on the scale s. 

Searching wider ranges of frequencies increases the value 
of /imax because the absolute number of gravita rs in w ider 
ranges of frequency increases. However, in (11.21 1 the 



gain from going to higher frequencies grows slowly, as 
the square root of the logarithm, because the gravitars 
spend less time at higher frequencies. 



C. A Natural Limit to the Result 

Because of the crucial assumption of a steady-state 
distribution in frequency there are obvious limits to this 
simple model. The time to reach a steady state in a given 
narrow frequency band [/, / -I- d/] must be at least of 
the same order of magnitude as the spin-down timescale 
(II.4 1, because otherwise no gravitar will have spun down 



to frequencies contained in the band. 

If all gravitars are born at the same high frequency the 
time to reach a steady state is exactly the spin-down age. 
If on the other hand there is a continuous distribution 
of initial frequencies, then reaching a steady state in a 
certain frequency band requires longer evolution times. 
Only then most of the gravitars in that band are ones 



that have spun down from higher frequencies. Over time 
this effect "washes out" any effects of the initial frequency 
distribution. 

There is a natural limit to the result due to the finite 
age of the Universe, since no gravitar can be older than 
the Universe itself. An even better limit would be the 
age of the Galaxy, or rather that of the galactic neutron 
star population. However, since the age of the Universe is 
known much more accurately than the age of the Galaxy, 
and since they differ only by a factor of order 2, we will 
use the age of Universe in all of our estimates below. 

The age of the Universe to can be calculated from Hub- 
ble 's constant Hq as io = f^o"^- Then the finiteness sets 
limits on the values of e and / for which the population 
has reached a steady state. We easily find from (n.4 1 that 



the population is in a steady state for gravitational-wave 
frequencies that satisfy 



s'f > ^HoP. 



(n.24) 



Fixing the ellipticity, one can calculate a frequency 

(n.25) 



/(.(.fHw^y.^H/'"-'^' 
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above which the population can be assumed to be close 
to a steady state. 

We would like to stress that a realistic population with 
a continuous distribution of initial frequencies has to have 
evolved over a time T « few x Tqw to be in steady-state. 
Thus the true value of / (e) is larger by a factor of a few, 
and falls into the frequency range of highest sensitivity in 
modern interferometric gravitational-wave detectors (be- 
tween 100 Hz and 300 Hz). 

The range of ellipticities for which the assumption of 
steady state breaks down is then given by 



Hf) ^5.8x 10" 



100 Hz 



(11.26) 



In general it is not valid to assume that the frequency 
distribution in our galaxy is in steady-state. 

We postpone further discussion of a uniform two- 
dimensional spatial distribution to Sec. |IVC| after pre- 
senting the setup of our numerical galactic model. 



III. NUMERICAL MODEL 

The second section of this paper gave a precise analytic 
formulation of Blandford's argument including the exten- 
sions and improvements of Ref . [? . To understand if this 
argument holds in a more realistic model of our galaxy, 
we set up a numerical simulation of the time evolution 
of a population of galactic gravitars. This follows [3], 
using a more recently published result [5] for the initial 
velocity distribution of neutron stars. To compute the 
spatial distribution, the equations of motion following 



TABLE L Mass and scale parameters for tire galactic potential 



Disk: 

Spheroid: 

Halo: 




= 8.07 X 101°Mq Od = 3.7kpc b^ 
= 1.12x10^°Mq as = 0kpc6s 
= 5.00 X 10^° Mq r„ = 6.0 kpc 


= 0.200 kpc 
= 0.277 kpc 



from the galactic potential given in Sec. |III A| are evolved 
over time. The assumed initial conditions for the differ- 
ential equations (i. e. initial positions and velocities of 
the gravitars) are described in Sees. IIIB and IIIC re- 



spectively. Sec. |III D] describes the adopted distributions 
for the initial spin period. The results of the simulations 
will be presented afterwards in Sees . |I V X| and [TV B [ 



A. galactic Potential and Equations of Motion 

The motion of galactic gravitars is governed by the 
galactic gravitational potential. The potential first given 
by Paczynski [6] is adopted. This potential describes our 
galaxy as axisymmetric with respect to the rotation axis. 
Thus, cylindrical coordinates p, z and ip are used, p 
denotes the distance to the galactic rotation axis, and z 
is the distance perpendicular to the disk. 

The adopted potential represents our galaxy as com- 
posed of three different mass components. The most 
massive is a nonuniform flat disk with a radial scale of 
3.7 kpc and a z-direction scale of 0.2 kpc. The compo- 
nent with the second highest mass is the halo, which 
is described by a density distribution Qh oc (r^ + r|) , 
where Th — 6 kpc is called the halo core radius. The cen- 
tral bulge of our galaxy is represented by a spheroidal 

_ 5 

lower mass component with a density g^ ex (r^ + &s) ^ i 
where 6s = 0.277 kpc. 

The corresponding potential therefore consists of three 
terms 



$ (p, z) = $s (r) + $D (p, z) + *H (r) . 



(IH.l) 



describing, respectively, the potential energy per unit 
mass of the spheroid, the disk and the halo in our galaxy. 
The first two components (i = S, d) are given by 



$, (p, z) = ~GM, 



P 



z2 + h] 



[111.2) 

For the potential of the halo r^ = p^ + z^ is substituted 
and the following spherical symmetric expression is used 



$H(r-) 



GMh 



In 1 



— arctan — 
r \ j'n 



(IH.3) 
The parameter values are shown in Table llj 

The axial symmetry of the galactic model leads to con- 
servation of the z-component of the angular momentum 



Lz- Thus, the effective potential is 



C. Initial Velocity 



$off(/0,2:) = <^{p,z) 



2p2 



The equations of motion that are evolved are 



P 



and z — —- 



dp 



dz 



(III.4) 



(III.5) 



The equation of motion for Lp is given by p^ip — L^. In 
the simulation this equation is not used because ^p is not 
evolved, but drawn from a uniform random distribution 

(p€E[0,27r). 



B. Initial Spatial Distribution 

The initial spatial distribution of gravitars is propor- 
tional to the density of massive progenitor stars of neu- 
tron stars. While there is quite good agreement about 
the initial distribution in z-direction, the initial distribu- 
tion along the radial direction is unknown. 

In the z-direction the initial position is drawn from a 
Laplacian distribution with scale factor zq = 0.075 kpc. 
The probability of a gravitar's birth in an interval 
[z, z -I- dz] is given by 



p^ (z) dz = -— cxp 

AZq 



2^0 



(III.6) 



We considered three different models for the initial 
radial distribution together with the given initial z- 
distribution. 

The simplest radial distribution [from Ref. [T] is an 
exponential fall off with a scale factor px — 3.2 kpc. The 
probability of a gravitar's birth in a distance interval 
[p, p + dp] is then 

Pi(p)dp=-expf-^')dp. (III.7) 

Pi V Pi/ 

Note however that this distribution leads to an extreme 
concentration of neutron stars towards the galactic cen- 
ter. These are not seen in pulsar surveys |7] . So following 
[B] a gamma distribution given by 



V2 (p) dp == a 



'Pl 



exp 



-^)dp, 

P2/ 



(III., 



is also considered, where gravitar formation in the disk 
is allowed for p ^ 25 kpc, and the constants are given by 
P2 = 4.78 kpc and Op = 1.0345. A third distribution 



P3 (p) dp : 



120p§ 



-^ dp, 

P3/ 



(III.9) 



with p3 = 1.25 kpc taken from [7^ is also considered. It 
is fitted to the radial distribution of Population I stars 
which are considered to be likely progenitors of neutron 
stars. Again, gravitar birth events are allowed for r ^ 
25 kpc. On average, only one gravitar out of 14 thousand 
is born with p > 25 kpc, so the normalization constant 
ao w 1. 



The galactic rotation determines the velocity of the su- 
pernova progenitors and therefore also that of the new- 
born neutron stars. 

From classical mechanics the rotational speed of a body 
on a circular orbit in the axisymmetric potential (III.ll 
is given by Vrot = \J pdp'^ (P: z). From the initial coordi- 
nate of a gravitar the corresponding Uj-ot on a tangential 
circular orbit is calculated (neglecting initial rotation ve- 
locities perpendicular to the galactic disk, because all 
gravitars are born with low initial values of z). Looking 
down on the Galaxy from positive z-values the Galaxy is 
chosen to rotate counterclockwise. 

Furthermore, it is assumed that gravitars are born in 
a supernova explosion that will kick the newborn star. 
The direction of that kick is assumed to be isotropic. 
The kick speed is drawn from a Maxwellian distribution 
with a mean velocity v = 430km/s. The probability for 
the kick speed to be in an interval [wkickj "^^kick + di^kick] is 



Pv (wkick) dv 



kick 



— -3^ exp 



4^kick 1 d^,. 



kick- 



(III.IO) 
It is assumed that this distribution from Ref. [5j may be 
used for gravitars as well as for pulsars. 



D. Initial Period Distribution 

Because of the conservation of angular momentum in 
the supernova event and the much smaller radius of the 
gravitar compared with its progenitor, newborn neutron 
stars will spin rapidly. 

We considered three models for the distribution of the 
initial rotation periods following the models given in Ref. 
[i\. The existence of three different models reflects our 
ignorance of the actual distribution of initial periods. 

The first model is described by a lognormal distribu- 
tion 



PPo (Po) = 



1 



V2TraP( 



exp 



1 



(in(Po)-ln(Po))^ 



(III.ll) 

where Pq is measured in seconds and where the values 
CT = 0.69 and Pq = 5 ms are taken from Ref. [9 . Gravitars 
with Pq < 0.5 ms are excluded. 

The second model is using the same probability distri- 
bution as the first model but every initial period Pq < 
10 ms is set to 10 ms exactly. In this way the possible 
presence of r-modes in young neutron stars is mimicked. 
These modes can dissipate rotational energy of the new- 
born neutron star and increase its initial period to about 
10 ms within 1 yr. 

The third model considered is a further extension of 
the second one. It includes the effects of matter fall-back 
after the supernova explosion. The increase of angular 
momentum by the accreting matter could counteract the 



r-mode induced deceleration. The resulting initial period 
will approach an intermediate value. The choice from [3] 
to draw the initial period from a uniform distribution 
between 2 ms and 15 ms is adopted 



E. Coding and Implementation 



and the equations of motion (III.5I at our disposal, it is 



a straightforward problem to find the spatial distribution 
of a population of gravitars at the present time. 

The code for the simulation is written in C. The equa- 



tions of motion (III. 5 1 are integrated via a Burlisch-Stoer 



method in combination with Stoermer's rule for the di- 
rect discretization of a system of second-order differential 
equations using routines described in [TU]. Over the in- 
tegration time the total energy is conserved to one part 
in 10^. 



For the derivation of the frequency distributions (11.15) 



the probability distributions are implemented by random 
number generators and functions from the GNU Scien- 
tific Library (GSL) [8] . The fractional frequency density 
(11.22 1 is obtained via a Monte-Carlo integration using 
2 X 10^^ random values of initial frequency and a uni- 
form distribution of ages. 

Depending on the model of the initial spatial distribu- 
tion the integration of 10^ neutron star trajectories over a 
time of 200Myrs takes 2.5mins to 13.3 mins on an AMD 
Opteron 185 processor. Most of the simulations were 
done on the Morgane cluster at the AEI in Potsdam. 



IV. RESULTS 



A. Frequency Distributions 



In the derivation of the generalized result (11.19 1 the as 



sumption of a steady-state frequency distribution (11.23) 



is (along with the two-dimensional uniform spatial dis- 
tribution) the key to the independence of the maximum 
amplitude on ellipticity and frequency. Let us therefore 
first have a look at the frequency distributions that result 
from a continuous distribution of initial frequencies and 
compare them with the corresponding density resulting 
from a single birth frequency (11.23 I. 

Fig. [l] shows the results of the Monte Carlo integra- 
tion for t = 13.6 Gyrs using a lognormal distribution 
of initial periods and a fixed ellipticity for each run in 
the frequency range / G [50 Hz, 2000 Hz]. The ellipticity 
varies over 3 orders of magnitude from 10^^ to 10^^. The 
dashed lines correspond to a scaling proportional to /^^, 
which results from a single birth frequency and shows a 
steady-state distribution as calculated in Eq. (11.23). 

For e = 10~^ the population is close to a steady state 
at the present time because the scaling is nearly propor- 
tional to f~^. For smaller ellipticities one can identify 
a kink in the density function at a frequency / as given 
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FIG. 1: The distribution gf{e,f) in frequency after i = 
13.6 Gyrs for varying ellipticity using a lognormal distribu- 
tion of initial periods in the range / G [50 Hz, 2000Hz]. The 
dashed lines correspond to a slope /^*, as Eq. (11.23 1 would 
predict for a single high birth frequency. The kink in the 



graphs for e < 10 ^ is at the frequency given by Eq. (11.25 1 



by (11.25). The kink is due to gravitars born at high fre- 



quencies that are too young to have spun down to lower 
frequencies. The frequency distribution for e < lO"'' is 
not in a steady state in the frequency range of highest 
sensitivity for modern interferometric detectors, which is 
between 100 Hz and 300 Hz. 

We also note that Qf (e, /) does not scale as e^'^ in all 
frequency bands. It only scales as e~^ at high frequencies. 

Let us now compare the fractional frequency densi- 
ties that result from different models of initial frequency 
distributions for a fixed ellipticity. Fig. |2] shows the dis- 
tribution in frequency space after an evolution time of 
t — 13.6 Gyrs for e — 10^^. Each graph corresponds to 
one of the models for the initial frequency distribution. 

If the evolution time is long compared with the spin- 
down time Tgw (H-4), most of the gravitars will have 



spun down to low frequencies, and the distribution will be 
dominated by those older sources. Yet, if the respective 
model has upper or lower limits on the birth frequen- 
cies one cannot expect the fractional densities to agree 
near these boundaries. However, in frequency bands of 
interest for modern interferometric detectors (100 Hz to 
300 Hz) the distributions only show minor differences be- 
tween the different models. 



B. Spatial Steady-State Distribution and 

Timescales 



Since the gravitars are born in a thin disk and receive 
an isotropic kick by the supernova, they tend to leave 
the disk after some Myrs. They either escape the galac- 
tic gravitational potential, or are bound to the Galaxy 
on some "orbit" . The numerical simulation is used to 
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FIG. 2: The distribution in frequency space in the range / € 
[50 Hz, 4000 Hz] after t = 13.6 Gyrs for e = 10"^ for three 
different models of the initial gravitar frequency distribution. 



find the timescale on which these processes wash out the 
imprint of the initial spatial distribution. 

Let us introduce the function M (r, t) , which is the 
number of gravitars in a ball of radius r around the posi- 
tion of the Sun that were formed a time t ago. The radial 
probability distribution g^ (r, t) as introduced in Sec. II B 
is related to M via the derivative with respect to r 



ft (r, t) dr 



1 



Nt, 



-drM{r,t)dr. 



(IV.l) 



To obtain a dynamical picture, Af (r, t) is computed in 
steps of 1 Myr from Myrs to 200 Myrs. For each of the 
201 values of integration time the trajectories of iVtot = 
10^ galactic gravitars are evolved using the numerical 



integration methods described in Sec. HIE and the radial 
distance distribution is derived from their final positions. 
The radial resolution is chosen as 2.5 pc for Okpc ^ r < 
12kpc and as 100 pc for 12kpc ^ r ^ 20kpc 

Fig. ISlshows the number of gravitars inside balls of ra- 
dius r around the position of the Sun, which were formed 
t = 200 Myrs ago, for the different initial radial distribu- 
tions. The saturation near r — 8.5 kpc is due to the high 
density of gravitars near the galactic center. 

The simulations show that the spatial distribution set- 
tles into a state of equilibrium for t « 200 Myrs. There 
is no significant difference between the distributions for 
t = 200 Myrs and i = 2 Gyrs. By the age of 200 Myrs the 
initial distribution is washed out; evolution over longer 
times no longer changes the radial distribution M (r) . 



C. Scaling Dimension of the Spatial Distribution 

The assumption of a two-dimensional and uniform spa- 
tial distribution of gravitars at the present time is crucial 
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FIG. 3: The distribution of gravitars as a function of radial 
distance from the Sun. The graphs show the number of grav- 
itars inside spheres of radius r plotted against r. The age of 
all sources is i = 200 Myrs. The scaling dimension of these 
graphs is shown in Fig. |4] The straight (dashed) lines show 
the corresponding slope for a uniform two- (three-) dimen- 
sional distribution. 



for Blandford's argument. The numerical simulation can 
test if these assumptions are valid or not. 

A useful concept is that of the scaling dimension. To 
obtain the scaling dimension, the function M (r) as intro- 
duced in the previous section is used. Assume a uniform 
distribution, and describe the number of sources inside 
each ball as a function of its radius r by a simple power 
law 



M (r) ex r^ 



(IV.2) 



D is called the scaling dimension of the distribution. 

Even an exactly two-dimensional spatial distribution of 
gravitars (e. g. the galactic disk) can effectively scale as 
a /^-dimensional object due to density gradients. Note, 
that the scaling dimension is a local quantity depending 
on the position of evaluation. 

To illustrate the scaling properties of the evolved galac- 
tic spatial distribution of gravitars, the local scaling di- 
mension is computed for r ^ 8 kpc and is shown in Fig. [4] 
From (IV.2 1 the scaling dimension D (r) can be derived 
via 



Dir) 



rdrM (r) 
AI{r) 



(IV.3) 



The differentiation is computed numerically using a cubic 
splining on the tabulated values. 

The scattering of points for r < 1 kpc is due to the 
small number of sources at short distances and resulting 
numerical noise. 

The scaling dimension at every distance to the Sun is 
greater than 2. For the first model of initial radial dis- 
tributions it even reaches values D ^ 5. For the second 
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FIG. 4: The local scaling dimension D (r) of the spatial dis- 
tribution of gravitars around the position of the Sun as cal- 
culated by (IV.3I. The corresponding radial distributions of 
sources are shown in Fig. [S] Note, that here the graphs are 
only shown for r ^ 8 kpc. 



model there is a slight increase in the scaling dimension 
towards the galactic center where the scaling dimension 
D Ki 3. For the last model the scaling dimension de- 
creases with larger radial distance but always D > 2. 

We conclude that the scaling dimension of the popula- 
tion of gravitars in the model of our galaxy is significantly 
larger than 2. More precisely averaging D over distances 
r ^ 2 kpc - where D is nearly constant and independent 
of the initial radial distribution - yields D « 2.95. 



V. THE STRONGEST CONTINUOUS 
GRAVITATIONAL- WAVE SIGNAL 

Using the frequency and spatial distributions obtained 
from our galactic simulation, it is straightforward to 
derive the maximum expected amplitude of continuous 
gravitational waves from gravitars. 



A. Numerical Method 

Let us first describe the numerical method for comput- 
ing the maximum expected amplitude of the gravitational 
waves using the distributions in space and frequency as 
presented in Sees . |I V !5] and [TV B | 

Eq. (11.18) is used to obtain the value of 
-^(/i.e/i./imax) for a given frequency band [/i,e/i] 
(In (e) = 1) and a trial value of /imax- To compare the 
results with Ref. |2] /imax is tuned via a bisection method 
within ±2.5% to the target value Aftar = 0.5 such that 
0.4875 < M (/i, e/i, fe^ax)^ 0-5125. 

the function M {r, t) giv- 



As described in Sec. 



IV B 



ing the number of gravitars with a solar radial distance 
less than r is tabulated. From these values g^ (r, t) = 



drM (r, t) /Ntot is numerically computed via a cubic 
splining method. 

The distribution in frequency space gf (e, /, t) is taken 
from Eq. (11.15 1 with distributions of initial periods as 
described in Sec. IIIIDI 

Given a frequency band a high value of /imax is cho- 
sen as trial value. Then the frequency is chosen fixed 
at the lower boundary of the band and the integration 
over h g [/j-ma x, feup] in ( II.18| is conducted by calculat- 
ing from ( |II.5[ ) the corresponding r (h) and inserting into 
the interpolated m (r, t) . For the upper limit /i^p of the 
integration the value that would obtained if the gravi- 
tar was at the closest possible distance is taken. Then 
no contribution to the integral is lost. After integrating 
over h the integration over frequency in the chosen band 
[/i, e/i] is done and weighted by gf (e, f,t). 

The result of these two integrations is then inte- 
grated over all times i € [0,t= 13.6 Gyrs] where the 
timestep between two evaluations is chosen as dt = 
min {iMyr, ^} to obtain a sufficiently fine timestep to 
track both spatial and frequency evolution. This inte- 
gration over time is weighted by a constant birthrate of 
h{t)^ {30 yrs)-\ 



B. Maximum Expected Amplitude 

Using the results of the simulation, one can see if the 
value of /imax (obtained by the method described in the 
previous section) differs from the one predicted by Bland- 
ford's result as extended and improved in Ref. [2]. 

To avoid boundary effects, from now only the first 
model for the distribution of initial frequencies is con- 
sidered. In all frequency bands, the other models always 
give a smaller maximum amplitude of gravitational waves 

"^niax ■ 

Fig. [5] shows the resulting value of /imax for different 
values of the ellipticity and all spatial distribution mod- 
els. The graphs are to be understood as follows: the 
maximum amplitude of gravitational waves is calculated 
as described in the previous section for every choice of ini- 
tial spatial distributions and for frequency bands [/, e/] . 
The graphs shows the value /imax obtained in such a band 
as a single point at (/, /imax)- 

From Fig. |5] it is obvious that the assumptions of 
Blandford's argument are not fulfilled for a realistic 
model of our galaxy. The graphs for different values 
of ellipticity do not line up, and each single graph is 
curved. Thus, the maximum amplitude of gravitational 
waves from galactic gravitars does depend on both the 
ellipticity and frequency. 

For highly deformed gravitars the graphs are nearly 
flat indicating a weak dependence on /, while for low 
values of e the previously-discussed kink appears. 

The choice of the initial radial distribution causes only 
small differences. The inter-model differences are usually 
of the order of 10%, and in the worst case are about 50% 
for a particular ellipticity and frequency. This is because 
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FIG. 5: The maximum strain amplitude of gravitational 
waves /imax from galactic gravitars in frequency bands [/, ef] 
in the range / G [50 Hz, 1000 Hz] . Each plotted point ftmax (/) 
is the value /imax calculated for a frequency band [/, ef] which 
is one natural-logarithmic-octave wide. The three curves show 
different initial spatial distribution models. For contrast, the 
dotted, straight line (independent of frequency and elliptic- 
ity e) shows the value of /imax from [2], which improved and 
extended Blandford's argument. 




gravitational wave frequency in Hz 

FIG. 6: The expected distance to the strongest gravitar, 
corresponding to Fig. [S] Because each frequency axis value 
of Fig. [5] refers to an entire logarithmic octave ln(/2//i) = 1 
range, in principle the plots should show a range of distances 
which is two logarithmic octave wide: ln(/i2//ii) = 2. To 
simplify the appearance of this plot, we show only a single 
curve at the central value. The expected range of distances 
ranges from a factor of 1/e below this plot to e above this 
plot. 



TABLE II; Maximum values for the amplitude /imax of grav- 
itational waves in dependence on the ellipticity e. 



10" 
10" 
10" 
10" 



1.6 X lO"^'* 

6.6 X lO"^'^ 
1.5 X 10"^^ 

2.7 X 10^^^ 



frequency band 

[250 Hz, 680 Hz] 

[550 Hz, 1500 Hz] 

[1000 Hz, 2800 Hz] 

[1000 Hz, 2800 Hz] 



the strongest gravitational waves are emitted by gravitars 
at very small distances [Tl]. As can be seen from Fig. 3] 
the distribution for radial distances r < 5kpc from the 
Sun is nearly the same independent of the initial radial 
distribution. 

Assuming the highest possible ellipticity e = 10~^ for 
the gravitars, the strongest signal has an amplitude of 
/imax « 1.6 X 10^24 j^ ^^le band [250 Hz, 680 Hz] which 
is improved (strengthened) by a factor of « 3 compared 
with the value given in ^. Note that this improvement 
factor would be smaller if larger values of e were consid- 
ered. In any case, the improvement factor is substantially 
larger at other frequencies. 

In the case of the more realistic value e = 10^^, 
the estimate is /imax ~ 6.6 x 10~^^ in the band 
[550 Hz, 1500 Hz] and is lower than the simple analytic 
estimate by a factor of 6. 

Table [n] lists the maximum values for the amplitude of 
the strongest gravitational waves for all adopted values 
of the ellipticity and the frequency band in which the 
maximum amplitude is expected. 



To illustrate the expected distance to the strongest 
gravitars we used the results of Fig. |5] together with Eq. 
(II. 5 1 to compute the distance to these sources as a func- 



tion of the gravitational- wave frequency interval and the 
ellipticity. The result is shown in Fig. [6] 

We conclude that the assumptions of Blandford's ar- 
gument do not hold in our galactic model. The expected 
gravitar spatial distribution is not a two-dimensional uni- 
form thin disk, and the expected gravitar frequency dis- 
tribution is not yet in a steady state for realistic values 
of neutron star ellipticity. 



C. Remarks on an Upper Limit 

In the previous section we obtained the expected max- 
imum amplitude of a gravitational-wave signal from a 
deformed neutron star spinning down purely by gravita- 
tional waves. Let us now address the question whether 
this value poses an upper limit on the gravitational- wave 
amplitude from objects that spin down partly by gravita- 
tional waves and partly by electromagnetic dipolar emis- 
sion. It must be stressed that recycled millisecond pulsars 
are not covered here, since spin-up is not considered. 

Ref. [2] gave a clever argument about why the expected 
maximum gravitational-wave strain from gravitars sets 
an upper limit for all neutron stars that have not gone 
through an accretion-powered spin-up phase. However 
that argument implicitly assumes only birth frequencies 
above the observed frequency band for the pulsars, and 
also assumes that the pulsar population frequency-space 
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distribution is in steady-state in the observed frequency 
band. Here, we do not make either of these assumptions. 

In this section we show that the results from Fig. [5] are 
a strict upper limit for frequencies / > 250 Hz. This is 
related to the distribution of initial frequencies which has 
its maximum at /o = 250 Hz and is decreasing monoton- 
ically for larger values of initial frequency. To maximize 
the number of sources in a frequency band at a given 
time in this regime, the slowest possible spin-down is re- 
quired. Any faster spin-down only would remove sources 
from this frequency band without adding more new ones 
from higher frequencies. Since the slowest possible spin- 
down without weakening the gravitational-wave signal is 
attained by switching off dipolar emission, the results 
from Fig. [5] are a strict upper limit on the gravitational- 
wave amplitude for / > 250 Hz. 

To formulate a rigorous proof of this claim, we char- 
acterize the generalized spin-down such a "mixed" neu- 
tron star will experience with a spin-down parameter y. 
The spin-down from electromagnetic dipolar emission in 
terms of the rotation frequency v is given by 



2tt^ B^R^ sin^ (a) 



7dipl^ 



(v.i) 



where Bp is the magnetic field strength at the neutron 
star's magnetic pole, a is the angle between the rotation 
axis and the magnetic field, and R is the radius of the 
neutron star. Rewriting Eq. (II.ll the spin-down from 
gravitational waves is 



512tt^G o . , 



(V.2) 



A neutron star emitting energy by both mechanisms at 
once will experience a total spin-down 



V = 7gwi'^ +7dipi'^- 



(V.3) 



This differential equation cannot (in contrast to (II.ll) 



be solved analytically for ly (t) , yet it can be integrated 
to give the time t (y, vq) in which a neutron star spins 
down from rotation frequency ^'o to v. 



t {v, I/q) 



r,,2 



I7dip| 



,2„2 



yln 



wl 



1 



yv' 



(V.4) 

where y :— 7gw/7dip is the general spin-down param- 
eter. Note, that y ^ corresponds to switching off 
gravitational wave emission, while y —^ oo leads to pure 
gravitational-wave spin-down. Taking these limits in 
(V.4 1, one easily recovers the equations for pure dipo- 



lar and gravitational- wave spin-down, respectively, which 
can be solved analytically for t/ (t) . 

Let us now turn to the derivation of the expected 
maximum gravitational wave amplitude from such neu- 
tron stars. The derivation given in Sec. |IIB| is 
straightforwardly modified to incorporate the general- 
ized spin-down. The generalization of the present-time 



gravitational- wave frequency distribution gj is obtained 
by writing Eq. (11.13 1 for fixed t and y as 



gj (/, t, y) d/ = g^, (/o (/, t, y)) ^M^:M. d/. 



(V.5) 



Note, that (V.4) cannot be solved analytically for 
/o (/i ^: y) = 2i/o (/, i, y). Care has to be taken to evalu- 
ate the partial derivative. Taking the total derivative of 
t{i',i'f)) — const with respect to v and application of the 
chain rule yields by a straightforward calculation 



df dv dv 



dt 
dvQ 



(V.6) 



Evaluation of this expression by use of (V.4 1 finally leads 
to 



df 
It is easy to see that y 



ducing Eq. (11.15 I 



fU^ + yfi) 
IH^ + yP)' 

-^ oo implies ^ 



(V.7) 



41, repro- 



However, there is no conceptual difference between a 
spin-down governed by e and one governed by y. One 
can write Eq. (11.181 giving the number of neutron stars 



with fixed spin-down parameter y in a frequency band 
[/ij/2] and gravitational- wave amplitude h ^ /imax as 

A'/(/i,/2, Vax,y) = 



h 



dtn{t) j dfgfif,t,y)J°° dh g, (r (h) , t) ^^ . 

(V.8) 

The fact that /imax for given y, /i, /2 and M is an 
upper limit on the gravitational-wave amplitude from a 
population of neutron stars, whose spin-down is governed 
by y, can be rephrased as follows: /imax is an upper limit, 
if M (/i, /2, /imax, y) is maximal as a function of y. If M 
is not maximal as a function of y, then a larger value 
of ft-max in the same frequency band could be found for 
a different value of y giving the same M. Thus, it is 
necessary to identify maxima of M (/i, /2, /imax, y) in y. 
These satisfy 

d 



dy 



M ifi,f2,h^iiK,y) = 0. 



(V.9) 



From Eq. (V.8) it is clear, that the only term affected by 
the derivative is gf {f,t,y). Applying the chain rule to 



(V.5 1 after inserting (|V.7[) one obtains 



dy 



3f 



dfo ' dy ■ fHA + yp) 



Qfo 



Jo [Jo J 



Pi'^ + ypy 

(V.IO) 
where the arguments of the functions are suppressed. 

Let us show that the last equation implies that the 
maximum expected amplitude as given in Sec. VB is a 
rigorous upper limit for frequencies /o for which 

dQfo 



dfo 



\foU,t,y) sSOV/,t. 



(v.ii) 



The second summand on the right-hand side of ( V.IO I is 
positive y f,t,y. The last factor in the first summand is 
always positive. It is also clear, tha t dyfp < V /, < for 
a fixed value of e [TS]. Therefore, if ( V.ll I holds true, it 



follows ^gj >0\/ f,t,y and ^Af > V y. Thus, if 



dy 



(V.ll I is fulfilled, the global maximum of M as a function 
of y is reached for y —^ oo, that is for pure gravitational 
wave spin-down. 

As mentioned earlier, for the model distribution of ini- 
tial frequencies adopted to produce Fig. |5] has its maxi- 
mum at /o — 250 Hz and is monotonically decreasing for 
larger values of initial frequencies. Therefore, Eq. (V.ll) 
is fulfilled for / ^ 250 Hz in Fig. |5] and the graphs shown 
there are a strict upper limit on the gravitational-wave 
amplitude from neutron stars spinning down by gravi- 
tational waves and electromagnetic dipolar emission for 
/>250Hz. 



VI. CONCLUSIONS 
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pie geometrical reasoning behind Blandford's argument 
is not valid. The numerical simulations provide an im- 
proved estimate of the expected maximum amplitude of 
gravitational waves from gravitars. We also showed that 
for frequencies / > 250 Hz this maximum amplitude is 
an upper limit for gravitational waves from neutron stars 
that spin down by gravitational waves and electromag- 
netic dipole emission. 

Although the expected maximum amplitude is lower 
by about 1 order of magnitude compared to the previous 
estimates, we would like to stress that in all of the models 
so far the influence of the Gould belt has been neglected. 
This young star- forming region (age ~ 40 Myrs) near the 
Sun is characterized by an abundance of massive O- and 
B-type stars enriching the solar neighborhood with young 
neutron stars. If there exists a population of gravitars 
born in the Gould belt, their gravitational-wave signals 
are more likely to be the first ones to be detected. 



We have used analytical arguments and the results of 
a numerical simulation to show that the assumptions of 
Blandford's argument do not hold in a realistic model of 
our galaxy. 

The assumptions (both in the original and in revised 
formulations of the argument) cannot be fulfilled for re- 
alistic values of ellipticity. The spatial scaling dimension 
D of an evolved neutron star distribution fulfills D > 2 
making a simple two-dimensional model invalid. The dis- 
tribution in frequency will not be in a steady state at the 
present time for realistic values of ellipticity. 

Because these two assumptions do not hold, the sim- 
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